;function_name must specify a function of form:
; name(neighbours) or
; name(neigbours, image[i,j])

;stupid circular kernel that simply tests each element in a square kernel for membership
function sjwkernc, image, radius, function_name, do_edges=do_edges, edge_policy=edge_policy
	sz =size(image, /dimensions)
	output = image
	output[*] = 0
	if (radius/2d eq radius/2l) then message, "Radius must be odd"
	if keyword_set(do_edges) then begin
		xd = 0
		yd = 0
	endif else begin
		xd = radius
		yd = xd
		sz[0]-=radius
		sz[1]-=radius
	endelse
	celladx = rebin(lindgen(radius * 2), radius * 2, radius * 2) - radius
	cellady = transpose(celladx)

	;simply check for each cell in the window if it falls in the radius (euclidean distance check)
	dist = euclidean_distances(transpose([[celladx[*]], [cellady[*]]]), reform([[0,0]],2,1))
	use = where(dist le radius)
	celladx = celladx[use]
	cellady = cellady[use]
	for j=yd, sz[0] - 1 do $
		for i=xd, sz[1]-1 do $
			output[i,j] = call_function(function_name, image[celladx[*] + i, cellady[*] + j], image[i,j])
	return, output
end

function sjwkern, image, wsize, function_name, do_edges=do_edges, edge_policy=edge_policy, circular=circular
	if keyword_set(circular) then return, sjwkernc(image, wsize, function_name, do_edges=do_edges, edge_policy=edge_policy)
	sz =size(image, /dimensions)
	output = image
	output[*] = 0
	if (wsize/2d eq wsize/2l) then message, "Window Size must be odd"
	if arg_present(do_edges) then begin
		xd = 0
		yd = 0
	endif else begin
		xd = wsize/2l
		yd = xd
	endelse
	celladx = rebin(lindgen(wsize), wsize, wsize) - wsize/2L
	cellady = transpose(celladx)
	for j=yd, sz[0] - 1 do $
		for i=xd, sz[1]-1 do $
			output[i,j] = call_function(function_name, image[celladx[*] + i, cellady[*] + j], image[i,j])
	return, output
end

